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ABSTRACT 

The 21 cm emission and absorption from gaseous halos around the first gen- 
eration of stars substantially depend on the Wouthuysen-Field (W-F) coupling, 
which relates the spin temperature with the kinetic temperature of hydrogen gas 
via the resonant scattering between Lya photons and neutral hydrogen. There- 
fore, the existence of Lya photons in the 21 cm region is essential. Although the 
center object generally is a strong source of Lya photons, the transfer of Lya 
photons in the 21 cm region is very inefficient, as the optical depth of Lya pho- 
tons is very large. Consequently, the Lya photons Uq from the source may not be 
able to transfer to the entire 21 cm region timely to provide the W-F coupling. 
This problem is especially important considering that the lifetime of first stars 
generally is short. We investigate this problem with the numerical solution of the 
integrodifferential equation, which describes the kinetics of Lya resonant photons 
in both physical and frequency spaces. We show that the photon transfer process 
in the physical space is actually coupled to that in the frequency space. First, 
the diffusion in the frequency space provides a shortcut for the diffusion in the 
physical space. It makes the mean time for the escape of resonant photon in op- 
tical depth r media roughly proportional to the optical depth r, not r 2 . Second 
and more importantly, the resonant scattering is effective in bouncing photons 
with frequency v ^ u back to u . This process can quickly restore z/ photons 
and establish the local Boltzmann distribution of the photon spectrum around 
uq. Therefore, the mechanism of "escape via shortcut" plus "bounce back" en- 
ables the W-F coupling to be properly realized in the 21 cm region around first 
stars. This mechanism also works for photons injected into the 21 cm region by 
redshift. 

Subject headings: cosmology: theory - intergalactic medium - radiation transfer 
- scattering 
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1. Introduction 

The resonant scattering of Lyo photons with neutral hydrogen atoms leads to a local 
Boltzmann distribution of the photon spectrum around the Lja frequency uq with color 
temperature equal to the kinetic temperature T of hydrogen gas. Consequently, the spin 
temperature T s of the hyperfine structure of neutral hydrogen will be coupled to the ki- 
netic temperature of hydrogen gas. This is the so-called Wouthuysen-Field (W-F) coupling 
(Wouthuysen 1952; Field 1958, 1959). The W-F coupling is crucial to estimate the redshifted 
21 cm signal from the halos of first generation of stars, because the deviation of the spin tem- 
perature from the temperature of cosmic microwave background (CMB) T cm b is considered 
to be mainly caused by the W-F coupling (e.g. Furlanetto et al. 2006). 

The lifetime of the first stars is short. The ionized and heated halos around the first 
luminous objects are strongly time-dependent. The time scale of the evolution of the ex- 
pected 21 cm emission/absorption regions can be as small as 10 5 years (Cen 2006; Liu et al. 
2007). Therefore, the 21 cm signal models based on the W-F coupling would be reasonable 
only if the time scale of the onset of the W-F coupling is less than that of the 21 cm region 
evolution. This time evolution has been studied very recently (Roy et al. 2009a, hereafter 
referred to as Paperl). It concludes that the local Boltzmann distribution can form within 
a time scale shorter than 10 5 years, but it would take 10 5 years or even longer to reach 
its saturation (time-independent) state. Therefore, it is legitimate to assume that the W-F 
coupling is taking place in the 21 cm emission/absorption regions, but the intensity of the 
photon flux is less than the time-independent solution. 

Paperl studied, however, only the case of spatial homogeneity and isotropy. It is equiv- 
alent to assume that the Lya photons are uniformly distributed in the entire 21 cm signal 
region. This assumption is not trivial. All the Lya photons in the 21 cm regions come 
from first stars, either from direct emission of Lja photons, or from the Hubble redshifted 
photons. On the other hand, the 21 cm regions are highly opaque for Lya photons. The 
W-F coupling may not be uniformly available in the 21 cm signal region, if the time scale of 
the transfer of Lja photons in the physical space is longer than that of the evolution of the 
21 cm region. 

The problem of Lya photon transfer in optical thick media is not new. It has been 
addressed as early as 1960s in relation to the escape of resonant photons from opaque clouds 
(Osterbrock 1962; Harrington 1973; Avery & House 1968; Adams 1972). In these references 
it is shown that the diffusion of the photon distribution in the frequency space caused by 
the resonant scattering will be helpful to speed up the spatial diffusion. In the past decade, 
there are also many works on the escape of Lya photons from high redshift objects with 
the mechanism of Hubble redshift (e.g. Miralda-Escude & Rees 1998; Loeb & Rybicki 1999; 
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Zheng & Miralda-Escude 2002; Haiman & Cen 2005; Tasitsiomi 2006; McQuinn et al. 2007). 
The Hubble redshifted Lja photons are also easy to take spatial transfer. Although both the 
resonant scattering and the Hubble redshift are useful to solve the problem of Lja photon 
transfer in the 21 cm region, they rely on the change of photon frequency from u to u ± Aia 
Therefore, the local Boltzmann frequency distribution will be disturbed, even if it initially is 
in the state of local Boltzmann frequency distribution. Thus, it is unclear whether the W-F 
coupling keeps to work timely and uniformly in the 21 cm region. 

In this context, a time-dependent solution of the kinetics of Lya photons in both physical 
and frequency spaces of the 21 cm region is necessary. This is the topic of the current paper. 
We will show that the resonant scattering of Lja photons is effective to solve the problem of 
spatial transfer in optical thick 21 cm region as well as to keep the W-F coupling working. 
On the other hand, for 21 cm regions of short lifetime objects, the cosmic expansion does 
not provide effective mechanism for the spatial transfer of Lja photons. 

Similar to Paperl, we use the numerical solution of the time-dependent integrodifferential 
equation of the radiative transfer with resonant scattering. The numerical solver is based 
on the weighted essentially non-oscillatory (WENO) scheme (Jiang & Shu 1996). WENO 
scheme is effective in solving Boltzmann equations (Carrillo et al. 2003, 2006) and radiative 
transfer (Qiu et al. 2006, 2007, 2008). The algorithm related to resonant scattering has also 
been given in Roy et al. (2009b). This numerical solver has successfully passed the tests of 
analytic solutions and conservation of photon number. Therefore, it is a good candidate for 
computing the current problem. 

This paper is organized as follows. Section 2 addresses the physical problems of Lja 
photon transfer in the 21 cm emission and absorption regions of luminous objects. Section 

3 is on the time scale of the transfer of resonant scattering in the frequency space. Section 

4 presents Lja photon transfer in the frequency and physical spaces. These results can be 
used for the 21 cm signal region (§5). Discussion and conclusion are given in Section 6. The 
details of the numerical implementation are given in the Appendix. 

2. Radiative transfer problem of the 21 cm region 

2.1. Basic properties of the 21 cm region 

The property of the ionized and heated regions around an individual luminous object 
is dependent on the luminosity (or mass), the spectrum of UV photon emission, and the 
time-evolution of the center object (Cen 2006; Liu et al. 2007). We will not work on a 
specific model, but consider only the common features. These halos generally consist of 
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three spheres. The most inner region of the halo is the highly ionized Stromgren sphere, 
or the HII region, in which the fraction of neutral hydrogen /hi = ^hi/^h is no more than 
10~ 5 , where tihi and hh are, respectively, the number densities of neutral hydrogen HI and 
total hydrogen. The temperature of the HII region is about 10 4 K. The physical radius of 
this sphere is in the range of a few to a few tens of kpc. The second region is the 21 cm 
emission shell, which is just outside the HII region. The physical size of this shell is similar 
to that of the HII sphere. The temperature of hydrogen gas in the emission shell is in the 
range 10 2 < T < 10 4 K due to the heating of UV photons. The fraction /hi is in the range of 
0.1 to 1. The third region is the 21 cm absorption shell, which is outside the 21 cm emission 
region. The physical size is on the order of 100 kpc. The temperature of this region is lower 
than T cm b = Tcmb(1 + z), where Tqmb is the temperature of CMB today. The time scale of 
the formation of the halos is about 10 6 years. The lifetime of the halos is about the same as 
the lifetime of the first stars. 

At the epoch of redshift z ~ 20, the reionization region consists of isolated patches 
around first sources. Most Lya photons in the 21 cm regions should come from the central 
source and the subsequent re-emission processes. If one can estimate the center object 
as a Lya emitter, the emission of Lya photons in number per unit time would be about 
dN-Lya/dt = 10 53 s _1 . The recombination of HII and electron in the Stromgren sphere is also 
a source of Lya photons. If the physical radius of the Stromgren sphere is ~ 10 kpc, the 
emission intensity is about dN^ ya /dt = 10 49 s -1 . 

Using the parameters of the concordance ACDM model, the optical depths of photons 
with Lya resonant frequency uq in the 21 cm region are 

where do is the cross section of the resonant scattering at the frequency u , Ris the physical 
size of the considered sphere, and r is the distance R in the units of the mean free path of uq 
photons at redshift z. Eq.flIT) shows that the optical depth of the 21 cm signal regions with 
/hi > 0.1 is r(u ) > 10 6 . 

It is also useful to define a dimensionless time rj = crimaot as 

✓ rp \ 1/2 / 2Q \ 3 /n 022\ 
t = ,/oHnot, = 6-7 x 10-/ HI ' { W ) (^j (j^) V (2) 

This is the time in the units of mean free flight-time. For a time scale t ~ 10 5 yrs, the scale 
of i] is ~ 10 7 . 
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2.2. Problems with the W-F coupling 

There are, at least, two radiative transfer problems with the W-F coupling in the 21 
cm regions: A.) How to provide enough Lya photons in such opaque medium? B.) Can the 
radiative transfer in the 21 cm region keep the W-F coupling to work well? 

If photons always keep the frequency u , their spatial diffusion can be described as a 
random walk process (Chandrasekar 1943). The mean number of scattering required for the 
diffusion over a range R is on the order of r 2 . Thus, the time for the diffusion over the 21 
cm range is rj ~ r 2 , corresponding to the time t = T 2 /cn^jaQ > 10 12 years. This diffusion 
mechanism, obviously, is useless for the 21 cm regions. 

The time scale of the diffusion would be substantially reduced if the frequency of the 
Lya photons can have a small shift from u$ to u = Vq ± Au, because the optical depth r at 
the frequency u ± Au is significantly less than t(u ). One can estimate the frequency shift 
Au by the condition of optical depth t(u ± Au) ~ 1, which is 

(f)(x,a) 
0(0, a) 

where <p(x, a) is the Voigt function for the resonant line u profile as (Hummer 1965) 



<P(x,a) = ^ B dy- rj- — ^ ( 4 ) 

where the dimensionless variable x is defined by x = iu — Vq)/ Aup and Avp = u vt/c is the 
Doppler broadening of hydrogen gas with thermal velocity vt = ^Jk^T /2m}±. Therefore, x 
measures the frequency deviation Au = \u — uq\ in the units of the Doppler broadening. The 
parameter a in eq.(j4]) is the ratio of the natural to the Doppler broadening. For Lya line, 
a = 2.35 x 10" 4 (T/10 4 )- 1/2 . In terms of x, the solution of eq.© is x > 3. 

An effective mechanism of the frequency shift is given by the diffusion of the photon 
distribution in the frequency space. Considering this mechanism, the number of scattering 
required for diffusion over a physical distance R is no longer equal to r 2 , but roughly on the 
order of r (Osterbrock 1962; Harrington 1973; Avery & House 1968; Adams 1972). Since 
r oc R, the time scale of the spatial diffusion over size R is comparable to R/c. Problem A 
would then be solved. 

Problem B still remains. If photons with the frequency u < uq — Au or u > uq + Au take 
a faster spatial diffusion than the u photons, how can we keep the W-F coupling to work? 
Without u photons, one cannot have a local Boltzmann distribution around u . Therefore, 
the photons with the frequency u < uq — Au or u > uq + Au should be brought back to the 
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frequency vq. We must study whether the frequency space diffusion mechanism can restore 
Vq photons from photons with the frequency v 7^ u . 

Cosmic expansion also leads to a deviation of photon frequency from z/ to a lower one 
Uq — Au, which speeds up the spatial transfer. However, we also need a mechanism to restore 
Lya photons from Hubble redshifted photons. Therefore, in terms of the W-F coupling in 
the 21 cm region, we must study the kinetics of Lya photons in both physical and frequency 
spaces. 



3. Radiative transfer of resonant photons in the frequency space 

3.1. Equations 

We first estimate the time scale needed for the frequency shift. We can use the radiative 
transfer equation of a homogeneous and isotropically expanding infinite medium consisting 
of neutral hydrogen. The equation of the mean intensity J in terms of the photon number 
is (Hummer & Rybicki 1992; Rybick & Dell'antonio 1994) 

^ = -<j>(x,a)J(x,7i) 

f dJ 

+ / lZ(x,x';a)J(x',r))dx' + r Y— — \- S(x, rj). (5) 

The parameter 7 _1 measures the number of scattering during a Hubble time, given by 

J \n M J V - 022 / V 20 / 



The re-distribution function TZ(x, x') gives the probability of a photon absorbed at the fre- 
quency x 1 , and re-emitted at the frequency x. It depends on the details of the scattering 
(Henyey 1941; Hummer 1962; Hummer, 1969). If we consider coherent scattering without 
recoil, the re-distribution function with the Voigt profile eq.fll]) is 

lZ(x, x'\ a) = (7) 
1 



7T 



3/2 



e- 2 



\x-x'\/2 



— 1 / ^min U \ —if ^max U 

tan — tan 



du 



where x m i n = mm(x,x') and x max = max(x, s'). In the case of a = 0, i.e. considering only 
the Doppler broadening, the re-distribution function is 

K{x,x') = ^e 26:t '' +f ' 2 erfc[max(|x + 6|,|x , + 6|)], (8) 
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where the parameter b = hv^jmvTC = 2.5 x 10 _4 (10 4 /T) 1//2 is due to the recoil of atoms. 
The re-distribution function of eq.® is normalized as lZ(x,x')dx' = (p(x, 0) = <f>g(x) = 
-j=e~ x2 . The numerical algorithm to solve eq.([5l) has been given in Roy et al. (2009a, 2009b). 



3.2. Time scales of the frequency shift 

On the right hand side of eq. (jSJ), the first term is the absorption at the resonant 
frequency x, the second term is the re-emission of photons with frequency x by scattering, 
and the third term describes the Hubble redshift of photons. We first solve eq. (jSJ) by 
dropping the terms of absorption and re-emission. The equation is 

_ =7 _ + S(l ,,). (9 ) 

Assuming the source is S — C<p s (x), where (f) s (x) is the normalized frequency profile of the 
source photons. The analytic solution of eq.([9]) is (Rybicki & Dell'antonio, 1994) 

J(x,rj) = J(x + 7r/,0) + C^ 1 / cj) s (x')dx' (10) 

J X 

where J(x, 0) is the initial flux. Define the mean frequency by 

_ fxJ(r],x)dx 
J J(r],x)dx 

If the initial flux is J{x, 0) = 0, one can show that for any profile 4> s (x) we have 

x( v ) = -777/2. (12) 

As expected, the speed of redshift dx/drj is a constant. For the Hubble expansion, a frequency 
shift x needs a time 77 ~ ^~ l x. Thus, in order to have the frequency shift x > 3, the time 
scale is 77 > 6 x 10 6 , corresponding to t > 4 x 10 4 years. This scale seems to be short 
enough compared with the lifetime of first stars. However, Hubble redshift is less effective 
comparing with the resonant scattering. This point can be seen in Figure dJ which gives 
the solution of eq.© with the re-distribution function eq.®. The source is taken to be 
S(x) = (pg{x) = (l/v / 7r)e _x2 . We also take the parameter 7 to be 10~ 6 [eq.([6])], and ignore 
recoil, i.e. 6 = 0. 

Figure [1] shows that the diffusion in the frequency space leads to a flat plateau with 
width \x\ < 3 when rj is as small as rj ~ 10 4 , or t ~ 10 2 years. This time scale is much less 
than that of the Hubble redshift. Therefore, the major mechanism to produce photons with 
frequency \x\ > 3 is given by the resonant scattering. 



- 8- 



10 6 

10" 
^10° 
10 ! 
10 1 
10° 

-15 -10 -5 5 

X 



Fig. 1. — Solutions J(rj,x) of eq. (j5J) with the re-distribution function eq.(IH]). Parameters 
b = 0, and 7 = KT 6 . 

3.3. Bounce back mechanism 

From Figure [U we can see that the profile of photons are almost symmetric with respect 
to x = (or to v = uq) until about the time rj = 10 6 . The redshift effect can only be seen 
from the curve of i] > 10 6 . That is, the resonant scattering impedes cosmic redshift. The 
impediment is due to the "bounce back" of resonant scattering. Regardless of whether 
the frequency of the absorbed photons is larger or smaller than u , the mean frequency of 
the re-emitted photons is always u . Therefore, the resonant scattering will bring back some 
redshifted photons to the frequency u . Thus the net effect of resonant scattering (absorption 
and re-emission) is, on average, to bounce redshifted photons back to the resonant frequency 
uo, and to restore the symmetry with respect to x — 0. This bounce back mechanism is the 
key of restoring u photons and the W-F coupling (see §4). 

To illuminate the bounce back effect, we calculate the mean frequency x(rj) for the 
solutions shown in Figure [U The result is plotted in Figure EJ The solution of pure cosmic 
redshift eq. (fT2|) is also shown in Figure [2] as the dashed curve. We can see from Figure [2] 
that the speed of photon frequency shift dx/drj when considering resonant scattering is in 
general less than that in the case of pure cosmic redshift. When the time rj is large, many 
photons have been redshifted out of the frequency range of 4> g {x). In this case, the "bounce 
back" effect ceases, and redshift speed is recovered to dx/dr] = —7/2. From Figure[2]one can 
see once again that the Hubble redshift can produce frequency shift from x = to |x| ~ 3 
only when rj > 10 7 , or t ~ 10 5 years. This time scale may not be short enough to match the 
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Fig. 2. — x vs. 7?7 for the solution shown in Figure [T] (solid) and the solution of eq. (fl2"l) 
(dashed). 

21 cm region with a short lifetime. When the time scale of resonant scattering is less than 
the time scale of Hubble expansion, the frequency shift of the Hubble expansion slows down 
significantly by the resonant scattering. 



4. Radiative transfer (RT) of Lya photons in the 21 cm regions 

4.1. RT equation in spherical halo 

Considering a photon source located at the central region of a uniformly distributed 
expanding medium, we can use the RT equation of the specific intensity I(r],r,x, //) as 
follows 



ar m (i - oi a/ 

drj ^ dr r dfi ^ dx 

—cj)(x)I+ / 1Z(x, x')I{rj, r, x', ii)dx' d(j! + S, 



(13) 



where fi = cos 9 is the direction relative to the radius vector r. The dimensionless coordinate 
r is rescaled from the physical coordinate 



rp \ 1/2 



2.1 x KTV, 



3 f-1 



20 

T+z 



0.022 



r pc. 



(14) 
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With these variables, the propagation of a signal with the speed of light will be described by 
the equation r = rj + const. It would still be reasonable to use the isotropic approximation 
of the re-distribution (Mihalas et al. 1976). 

When the optical depth is large, the Eddington approximation would be proper. It is 

I(r], r, x, 11) ~ J(r), r, x) + 3/xF(r/, r, x) (15) 

where J(r/, r, x) = | j\ I(rj, r, x, fi)dfi is the angularly averaged specific intensity and F(rj, r, x) 
| pil(rj, r, x, fx)dfx is the flux. Defining j = r 2 J and / = r 2 F, Eq. fflBl yields the equations 
of j and / as 

% + ^ = -<t>{x)3 + [ Hi?, x')jdx' + 7^ + r 2 S, (16) 
Of] or J ox 

¥ + = -ms- (17) 

or) 3 or 3 r 

The numerical algorithm for eqs.( fT6l) and ( TTTI) is given in the Appendix. 



Since photons with frequency shifted away from vq would be optical thin, the equations 
( TT6|) and (fT7j) will no longer be a good approximation when the frequency of photons is 
shifted to the optical thin case. From Figures [1] and [2j one can see that within 77 < 10 6 , most 
photons are still trapped in the frequency \x\ < 3, for which the optical depth is larger than 
1. Therefore, the Eddington approximation would be proper, at least, until r\ is as large as 
about 10 6 . 



4.2. Diffusion in the physical space 

We first solve the equations (Tl6|) and (TT7T) by dropping all terms on the transfer in 
frequency. It has been shown that the source term can be replaced by a boundary condition 
of / and j at r = r (Qiu et al. 2006). The equations then become 

33 + ¥■ = 0, (18) 



% + ^-U= f (») 



drj dr 

dr] ' 3 dr 3 r 
We take the boundary condition at r = 1 to be 

j(77,l) = 3£ , f(v, 1) = So = 10 6 . (20) 
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Since eqs. ()18l) and f|T9|) are linear, the parameter So is not important if we are only interested 
in the shape of j and / as function of rj and r. The initial condition is taken to be 

j(0,r) = /(0,r) = 0. (21) 

The solution of f(rj,r) is presented in Figured The dashed line is / = S = 10 6 , which is 
the time-independently exact solution of eq. (IT5|) . It is a horizontal straight line because the 
flux / = r 2 F is r- independent and satisfies the conservation of the photon number. Figure [3] 
shows that the spatial size r of f(rj, r) roughly satisfies r oc y/rj. Therefore, without resonant 
scattering, the diffusion basically is a random walk process as discussed in §2.2. 
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Fig. 3. — The solution f(r),r) of eqs. (TT8j) and (fl9]l . The straightline is / = So — 10 6 . 

We now turn to the solution of eqs. ffTo 7 ]) and (jT7j) with resonant scattering and Hubble 
redshift. The relevant parameters are given by b = and 7 = 10~ 5 . We use the boundary 
condition at r = 

j( V ,0,x) = 0, f( V ,0,x) = S <f> a (x) (22) 

where the frequency profile is taken to be the Gaussian profile, i.e. <p s (x) = <j> g (x). The 
parameter So is still taken to be 10 6 . The initial condition is similar to eq. (l2~T|) . i.e. j(0, r, x) = 
f(0,r,x) = 0. The solutions of f(r],r,x) are plotted in Figure HI 

All the solutions of Figure H] show two remarkable peaks at x — ±(2 — 3) for all radius 
r. The amplitude of f(j],r,x) at the peaks is higher than that at x = by a factor of 10 
to 10 2 . It shows that the flux is dominated by photons with frequency x ~ ±(2 — 3). That 
is, the spatial transfer is carried out by photons of x ~ ±(2 — 3). The amplitude of the flux 
at the saturated peaks is basically r- independent, and the saturated values of / f(rj,r,x)dx 
are r-independent. This is consistent with the conservation of the photon number. 



Fig. 4. — The flux f(rj, r, x) of the solutions of eqs. fll6|) and f[T7|) at r = 50 (left), 100 (middle) 
and 500 (right). The frequency profile of the source is <p g {x) = (l/y / 7r)e~ a:2 . 




Fig. 5. — The flux f(r], r, x) of the solutions of eqs. (lT6i) and ( TT71) at r = 50 (left), 100 (middle) 
and 500 (right). The frequency profile of the source is s (x) = (1/ y / ir/2)e~ 2x2 . 



More importantly, the amplitude of the peaks approaches its saturation at about r] ~ 200 
for r = 50, 77 ~ 400 for r = 100, and 77 ~ 2, 000 for r = 500. That is, the time 77 needed for 
the spatial transfer over size r is roughly proportional to r. This is very different from the 
random walk relation i] oc r 2 , or the results shown in Figure El 

It should be emphasized that the photons at the two peaks are not only those with 
x ~ ±(2 — 3) directly from the source Sq(/) s (x), but also include the photons with frequency- 
shift from x ~ to x ~ ±(2 — 3). This point can be shown by a source of the profile with 
smaller width, say, <fi s (x) = (1/ y / n/2)e~ 2x2 . In Figure El we plot the results. For the source 
of (f) s (x) = (1/ a/7t /2)e~ 2x2 , the number of photons with x ~ ±(2 — 3) is much less than that 
of the source (l/^)e~ x2 . We can see that all the features in Figure [5] are the same as those 
in Figure HI The two peaks are still located at x ~ ±(2 — 3). The time scale for approaching 
saturation in Figure [5] is also the same as that in Figure HI 

Therefore, photons at the peaks of x ~ ±(2 — 3) should come from the frequency-shift 
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from x = to x ~ ±(2 — 3) by resonant scattering. That is, resonant scattering provides 
a shortcut of the spatial transfer: first to shift x = photons to x ~ ±(2 — 3) photons, 
and then to speed up spatial transfer. This mechanism is the same as that for the escape of 
resonant photons from opaque clouds (Osterbrock 1962; Harrington 1973; Avery & House 
1968; Adams 1972). It allows the Lya photons to be able to transfer over to the 21 cm 
regions in time rj proportional to its size r. 



4.3. Resonant photon restoration and W-F coupling onset 

We have mentioned in §4.2 that the flux f(rj,r,x) is lacking u Q (or x — 0) photons 
even when / approaches its saturation. This is, obviously, not good for the W-F coupling. 
However, we found that the situation may not be so if we consider the solution of the mean 
specific intensity j. Figure [6] presents the solution j(r],r,x) of equations (fTBT) and (jTTI) at 
r = 50, 100 and 500. The parameters 6 = and 7 = 10~ 5 are the same as the solutions in 
Figure HI The results are given in Figure EJ 




Fig. 6. — The mean intensity j(t],r,x) of eqs.( fT6|) and ( JT7l) at r = 50 (left), 100 (middle) 
and 500 (right). 

When the time 77 is small, j(r), r, x) is similar to f(t], r, x), having a valley around x = 
and two peaks at x ~ ±(2 — 3). However, different from the flux f(r], r, x), the amplitude of 
j(r],r,x) around x = is quickly increasing. In the saturated state it is about the same as 
the peaks. That is, although the flux always has a valley around x = 0, the vq photons are 
quickly restored in the mean intensity. The time scale of approaching its saturated state is 
also proportional to r, not r 2 . Therefore, the restoration of resonant photons is due to the 
resonant scattering "bounce back" (§3.3), which pushes photons with frequency x ~ ±(2 — 3) 
back to x ~ 0. 

The shape of j(r], r, x) around x = is a flat plateau, which is similar to that in Figure 
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[TJ As have been shown in Paperl, the flat plateau of j(rj, r, x) at b = will become the local 
Boltzmann distribution if the recoil is considered. We can expect that the flat plateau of 
Figure M will also show a local Boltzmann distribution if b ^ 0. We calculate the solutions 
j(rj,r,x) and f(r),r,x) of eqs. (fTBT) and ffTTl) at r = 500 with the same parameters as those 
in Figure 6, but with b = 0.3. We use a large b, because it is easier to see the slope of the 
local Boltzmann distribution. The results are given in Figure [71 

Figure [7J clearly shows a local Boltzmann distribution within the range |x| < 2 as 

r, x) ~ j(r], r, 0)e~ 2bx = j(rj, r, o) e -H^o)/k B T _ ^ 

We find the slope to be In j(r], 0) — In 1) = 0.59, which is well consistent with 2b = 0.6. 
Figure [7J shows that at r = 500, the onset of the W-F coupling can occur as early as r\ = 900, 
but the amplitude of the local Boltzmann distribution at that time is much lower than its 
saturated value by a factor of 10 2 . The amplitude of the local Boltzmann distribution is 
substantially increasing with time. Therefore, the "bounce back" mechanism keeps the W-F 
coupling to work with a timescale t] larger than a few hundreds, i.e. a few hundred collisions. 
This result is the same as that in Paperl. 
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Fig. 7. — The mean intensity j(r],r,x) (left) and flux f(j],r,x) (right) of the solution of 
eqs. (1161) and ffTTl) at r = 500. The recoil parameter is taken to be b = 0.3. 

The solution of the flux f(r], r, x) in Figure [7J is not very different from that in Figure HI 
The only difference between the two figures is that the former is asymmetric with respect to 
x = 0, i.e. the peak at x < is stronger than that of x > 0, while the latter is symmetric. 
This is simply due to the recoil of b ^ leading more photons to move to x < 0. Neither 
flat plateau nor local Boltzmann distribution is shown in the flux f(i], r, x). There is always 
a valley around z/ even when f(r/,r,x) is in its saturated state. It once again indicates that 
the flux is dominated by photons of x ~ ±(2 — 3). 
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Fig. 8. — The mean intensity j(r),r,x) (left) and flux f(r),r,x) (right) of the solution of 
eqs. (jT6l) and f|T7|) with the Voigt profile eq.(@| and re-distribution function eq.([7|). The 
parameters are taken to be r = 50 and a = 10~ 2 . 

The two components / and j of the Eddington approximation are effective in revealing 
the functions of the resonant scattering. The flux f(j],r,x) describes the photons in transit. 
It shows that the spatial transfer within opaque media is mainly via photons with frequency 
shifted to x ~ ±(2 — 3), which are easy for escaping. The mean intensity j(rj, r, x) describes 
the restoration of x = photons and the onset of the W-F coupling. Therefore, f(j],r,x) 
shows a deep valley around x = (or i/q), while j(rj,r,x) shows a plateau. 

In Figure [HI we present the the solutions of mean intensity j(j],r,x) and flux f(j],r,x) 
given by eqs. (|T6|) and (|T7|) with the Voigt profile eq.(jl]) and re-distribution function eq.©. 
As expected, Figure M has the Lorentz wings. However, in the center part \x\ < 3, Figure M 
shows the same features as Figures H] and O That is, the mechanism of "escape via shortcut" 
plus "bounce back", which mainly relies on photons with \x\ < 3, still works well. The 
Lorentz wing only leads to long tails in the profiles of j and / in the range \x\ > 3, and the 
wings have very low amplitudes. 

4.4. Effect of injected photons 

Hubble redshift is another mechanism to produce photons with frequency ~ z/ , which 
also have the problems of the spatial transfer and the W-F coupling. Since the 21 cm region 
basically is optical thin for photons with x > 3, we first model the cosmic redshift by a 
photon source with the profile (p s (x) = 4> g (x — 3). The solutions of f(r],r,x) and j(r],r,x) 
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Fig. 9. — The mean intensity j(rj,r,x) and flux f(j],r,x) of eqs.( fT6l) and (JTTJ) at r = 50 
(left), 100 (middle) and 500 (right), parameter b = 0.3. The frequency profile of the source 
is 4> s (x) = 4> g (x - 3). 

with the same parameters as in Figure [7] are shown in Figure [H 

Although the injected photon has a peak at x = 3, the flux f(r),r,x) in Figure [9] still 
shows two peaks at x ~ ±(2 — 3). The peak at x — (2 — 3) is higher than that at — x ~ (2 — 3). 
It is simply because the photon source is at x = 3. The peak at x < is also much higher 
than f(r),r,x) at x — 0. That is, even though the source photon is at x = 3, the spatial 
transfer is still dominated by photons of both x ~ (2 — 3) and x ~ —(2 — 3). The mean 
intensity j(j],r,x) shows once again a perfect local Boltzmann distribution with slope 2b in 
the range —2 < x < 2 (eq. (|23p ). Any photons injected by the redshift into the 21 cm region 
will quickly join the W-F coupling by the bounce back mechanism. 

We now model the source with a continuous frequency spectrum: <p s (x) = Sq within 
\x\ < 10, and <f) s (x) = at \x\ > 10. The solutions of the mean intensity j(r],r,x) and 
flux f(r],r,x) of eqs.( [T6l) and ( 1T71) are shown in Figure fTUl in which the other parameters 
are the same as those in FigureEl We can see that in the center part \x\ < 3, j(r],r,x) 
and f(rj,r,x) have the same features as in Figure This result is expected, as the center 



Fig. 10. — The mean intensity j(r), r, x) and flux f(r], r, x) of eqs. fTEI) and ( TTT1) with the Voigt 
profile eq.(jl]) and re-distribution function eq.([7]). The frequency profile of the source is flat, 
i.e. 4> s (x) is equal to S within \x\ < 10, and to zero at \x\ > 10. 

part has contributions from the photons redshifted to \x\ < 3 from x > 3. Figure [9] shows 
that the source of <fr g (x — 3) does not change the mechanism of "escape via shortcut" plus 
"bounce back", and therefore, the source with a continuous frequency spectrum will keep 
these features as well. 



4.5. Effect of the Hubble redshift 

Although the cosmic expansion is considered in all the above-mentioned solutions, the 
effect of cosmic expansion seems to be negligible. It is because the optical depth is large and 
the parameter 7 is very small. The number of resonant scattering within the Hubble time is 
very large. The cosmic expansion is too small in one free flight time. The "bounce back" is 
dominant. 

When /hi is smaller, optical depth of the z/ photons is smaller, and 7 is larger, the 
effect of the Hubble redshift would appear. Figure [TT1 presents the solution f(rj,r,x) and 
j(r),r,x) of eqs. ffT6l) and ffT7|) with the same parameters as those in Figures @] and [6], except 
that the parameter 7 = 10~ 3 is large. Both j(r], r, x) and f(r], r, x) show the same features as 
in Figures H] and The Hubble redshift makes the profile to be asymmetric with respect to 
x = 0. The red wing is stronger than the blue wing. j(rj, r, x) still shows a flat plateau in the 
range —2 < x < 2. Therefore, the W-F coupling will work when 7 = 10~ 3 , corresponding to 
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Fig. 11. — The solutions j(r),r,x) and f(r),r,x) of eqs. flTE]) and (JTTJ) at r = 50 (left), 100 
(middle) and 500 (right). 7 = 10~ 3 . The parameters are the same as those in Figures H] and 
IHl except that 7 = 10~ 3 . 

/hi ^ 10- 3 . 

5. An example of time-dependent effect of W-F coupling on 21 cm signals 

The time evolution of W-F coupling shown in §4 should be considered in calculating the 
21 cm emission and absorption from the halo around the first sources. There are many self- 
consistent models with different parameters, such as the UV photon luminosity, frequency 
spectrum, the index of power law spectrum, the temperature of black-body spectrum, etc. 
(e.g. Cen 2006, Liu et al. 2007, Chen & Miralda-Escude 2008). But no one of these models 
considered the time-dependence of the W-F coupling. 

To show the importance of the time-evolution of the W-F coupling, we re-calculate the 
brightness temperature 5Tb of 21 cm signals of one of the models developed in Liu et al. 
(2007) where source intensity is E = 7.25 x lO^ergs -1 as in their Figs. 3 and 4, which is 
self-consistent in terms of heating and cooling of gas. The results are presented in Figure 12. 
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Fig. 12. — The 21 cm brightness temperature 5T b as a function of time at three comoving 
distances r = 3, 6, 9 Mpc, respectively, for the UV heating model in Liu et al (2007). The 
dashed lines are for the solutions with time-independent W-F coupling, while solid lines show 
those in which the time dependence of the W-F coupling is considered. 



Figure 12 shows the time-dependence of the 21 cm brightness temperature at three 
shells with radius 3, 6, and 9 h _1 Mpc. For each radius, there are two curves, one does 
not consider the time-dependence of the W-F coupling, and one does. It shows that the 
brightness temperatures are affected by the time evolution of the W-F coupling significantly. 
The time dependence of the W-F coupling will make the 21 cm signals weaker, especially for 
absorption features by a factor 2-5. This is because absorption features are formed just after 
the light front passes a particular location, at a time when the intensity of Lya photons in 
the Boltzmann distribution is far lower than their saturated values (Fig. 6). The absorption 
areas of the 21 cm signal would suffer more from insufficient Lya photons because they are 
too close to the light front. On the other hand, the decrease in ST b is negligible for emission 
areas. 
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6. Discussion and conclusions 

The kinetics of Lya resonant photons in the HI media with high optical depth r can 
basically be described as diffusion in both the physical space and the frequency space. If 
Lya photons do not join the diffusion in the frequency space, the transfer of Lya photons 
in the physical space is very inefficient, as the number of scattering needed for escape is 
proportional to r 2 . The resonant scattering of Lya photons and neutral hydrogen makes the 
diffusion processes in the physical space coupled to the diffusion processes in the frequency 
space. First, the diffusion in the frequency space provides a shortcut for the diffusion in 
the physical space. It makes the mean number of scattering for escape to be approximately 
proportional to r. Second, the bounce back of resonant scattering provides a mechanism 
of quickly restoring vq photons from x ~ ±(2 — 3) photons. Finally the W-F coupling is 
realized simultaneously with the restoration of the x = photons. 

The mechanism of "escape via shortcut" plus "bounce back" is mainly carried out by 
the photons with frequency x ~ ±(2 — 3). In a 21 cm emission region of physical size R ~ 10 
kpc and /hi > 0.1, the optical depth of x ~ ±(2 — 3) photons is still larger than 1. Therefore, 
it is reasonable to use Eddington approximation. On the other hand, the optical depth of 
x ~ ±(2 — 3) photons is much less than that of the x = photons. The x ~ ±(2 — 3) 
photons can transfer and enter the 21 cm emission region in a time scale less than 10 5 years. 
Therefore, the mechanism of "escape via shortcut" plus "bounce back" is able to timely 
support the W-F coupling of the 21 cm emission shell with Lya photons from the center 
objects. 

The time dependence of the W-F coupling would make the 21 cm signals weaker than 
the predication given by models which do not consider this time-evolution. Especially at the 
early stage of the formation of the 21 cm signal regions, the intensity of the local Boltzmann 
distribution is still very low, and therefore, one cannot assume that the spin temperature of 
21 cm is locked to the kinetic temperature of gas. It may yield a low brightness temperature 
of the 21 cm signals. 

Although the mechanism of "escape via shortcut" plus "bounce back" helps Lya dif- 
fusion, it does not mean that this mechanism will reduce the Gunn-Peterson optical depth 
of the Lya photons. On the contrary, the resonant scattering will lead to a slight increase 
of the optical depth of the Lya in the 21 cm region, as the resonant scattering impedes the 
cosmic redshift (Figure [2]). Consequently, there should be no observable redshifted optical 
signal with the frequency vo/(l + z) to be spatially correlated with the (1 + z) redshifted 21 
cm signal. 

The evolution of photons described by eqs. tTTfij) and (TTTT) conserves photon numbers. The 
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number of Lyct photons is basically conserved if one can ignore the Lyct photon destruction 
processes, such as the two-photon process (Spitzer & Greenstein 1951, Osterbrock, 1962). 
Thus, subsequent evolution of the Lya photons in the 21 cm region is to diffuse to a large 
sphere around the first stars. At the same time, the Lya photons will be redshifted. When 
the redshift is large enough, their Gunn-Peterson optical depth will be small, and finally 
these photons will escape from the halo (e.g. Miralda-Escude & Rees 1998; Loeb & Rybicki 
1999; Zheng & Miralda-Escude 2002; Haiman & Cen 2005). The escaping sphere should be 
larger than the size of the 21 cm region. Therefore, redshifted Lya optical signal with low 
surface brightness may come from a big halo around the 21 cm emission region. 

This work is supported in part by the US NSF under the grants AST-0506734 and 
AST-0507340 and by ARO grant W911NF-08-1-0520. WX is grateful for the hospitality of 
National Astronomical Observatories of China, where part of this work was done. 



A. Numerical algorithm 

To solve equations (fl6l and ( fTTj) as a system, our computational domain is (r, x) G 
[0, r max ] x [xi e ft, bright], where r max , xi e f t and x rig ht are chosen such that the solution vanishes 
to zero outside the boundaries. In the following, we describe numerical techniques involved in 
our algorithm, including approximations to the spatial derivatives, integrals in the frequency 
domain, numerical boundary condition and time evolution. 



A.l. The WENO algorithm: approximations to the spatial derivatives 

The spatial derivative terms in equations (IT51) and (fT7j) are approximated by a fifth 
order finite difference WENO scheme. 

We first give the WENO reconstruction procedure in approximating |£, 

dj(rj n ,ri,Xj) 1 



dx Ax 



Oj+1/2 - V1/2), (Al) 



with fixed rj = r] n and r = r*. The numerical flux /ij+1/2 is obtained by the fifth order WENO 
approximation in an upwind fashion, because the wind direction is fixed (negative). Denote 

hj = j( V n , n, Xj ), j = -2, -1, • • • , N + 3 (A2) 

with fixed n and i. The numerical flux from the WENO procedure is obtained by 

hj+1/2 = wihfl 1/2 + + wshfli/2> ( A3 ) 
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where hv^^ are the three third order fluxes on three different stencils given by 



^•+1/2 - 3^' + g^' +1 ~~ g^'+ 2 ' 



1,5, 1 



and the nonlinear weights uo m are given by, 



Urn . 7« , A/1 x 

«m=— , ^=(7TaF' (A) 



i=l 

where e is a parameter to avoid the denominator to become zero and is taken as e = 1CT 8 . 
The linear weights ji are given by 

71 =10' 72 =5' 73 =10' (A5) 
and the smoothness indicators (3\ are given by, 

13 1 

A = J^(h j „ 1 -2h j + h j+1 ) 2 + -(Vi-4/i J + 3/ lj+ i) 2 , 

13 1 

A = ^ ( h j- 2h j+i + h j+2) 2 + ^(hj-h j+2 ) 2 , 

13 1 

A = —(h j+1 -2h j+2 + h j+3 ) 2 + -(3h j+1 -4:h j+ 2 + h j+3 ) 2 . 

To approximate the r-derivatives in the system of equations (fl6l) and (TT7|) . we need to 
perform the WENO procedure based on a characteristic decomposition. We write the left 
hand side of equations (fT6|) and (ITT)) as 

u t + Au r (A6) 

where u = (j, f) T and 

VI o 

is a constant matrix. To perform the characteristic decomposition, we first compute the 
eigenvalues, the right eigenvectors, and the left eigenvectors of A and denote them by, A, 
R and R~ l . We then project u to the local characteristic fields v with v = i? _1 u. Now 
Uj + Au r of the original system is decoupled as two independent equations as v t + Av r . 
We approximate the derivative v r component by component, each with the correct upwind 
direction, with the WENO reconstruction procedure similar to the procedure described above 
for In the end, we transform v r back to the physical space by u r = Rv r . We refer the 
readers to Cockburn et al. 1998 for more implementation details. 
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A.2. Normalization of the re-distribution term 

We apply the rectangular rule to evaluate J R(x,x')jdx. The rectangular rule is known 
to have spectral accuracy for smooth integrated function R(x, x')j(x') with compact support. 
In equations (fT6l) and ffTTl) . it is known that 

J R(x,x')dx' = (f)(x), (A7) 

which is crucial for photon conservation. However, 

R(x, Xj)Ax = (f){x) 

3 

may not be true in general. To numerically preserve the photon conservation, we first 
compute 

4>{x) 

ratio(x,Ax) = -r— , (A8) 

^jR^XjXjjAx 

then we normalize the collision term by approximating J R(x,x')j(t,r,x')dx' with 

ratio(x, Ax) R(x, Xj)j(r), r, Xj)Ax (A9) 

3 

with fixed r\ and r. 



A. 3. Implementation of the boundary condition 

The source term given in the equation f|T6l) is implemented as a boundary condition on 
f(r],r = r , a;). 

ffar = r Q ,x) = s (f) s (x) 

For the intensity j, a reflective boundary condition is used at r = r . At the boundary 
of r = r max , x = Xi e f t and x = x rig ht, we use zero boundary conditions for both j and /, 
because of the way we choose r max , X\ e ft and x r i g ht- 



A. 4. Time evolution 



The time derivatives and are approximated by the third-order TVD Runge Kutta 
time discretization (Shu & Osher, 1988). For systems of ODEs u t = L(u), the third order 
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Runge-Kutta method is 

= u n + AtL{u n ,f' 



u 



(2) 



3 „ 1 



-u n + + AtL{u {1 \t n + At)), 



4 4 



M n+1 = \u n + |(w (2) + AtL{u ( - 2) ,t n + ^At)). 



A. 5. Test with the conservation of the photon number 

From eq. (jlfip we have 

^Jjdx + -^-J fdx = (A10) 
Therefore, for time-independent solution we have J f(r,x)dx = const. It yields 

f(r,x)dx= J f(0,x)dx = So (All) 

That is, the flux at saturated states is r-independent. This is the conservation of the number 
of photons. It can be used to test the algorithm. 
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